//*********************************************
//  This is Geant4 Template
//                                  author:Qian
//

#include "MyG4BasedAnalysis.hh"
#include "g4root.hh"
#include "Verbose.hh"
#include <cmath>
#include "G4AntiNeutrinoE.hh"
#include "G4AntiNeutrinoMu.hh"
#include "G4AntiNeutrinoTau.hh"
#include <cstdlib>
#include "G4ProcessType.hh"
#include <ctime>
//fNotDefined, fTransportation, fElectromagnetic, fOptical, fHadronic, fPhotolepton_hadron,
//fDecay, fGeneral, fParameterisation, fUserDefined, fParallel, fPhonon, fUCN
#include "G4DecayProcessType.hh"
//DECAY, DECAY_WithSpin, DECAY_PionMakeSpin, DECAY_Radioactive, DECAY_Unknown, DECAY_External
#include "G4HadronicProcessType.hh"
//fHadronElastic, fHadronInelastic, fCapture, fFission, fHadronAtRest, fLeptonAtRest, fChargeExchange, fRadioactiveDecay
#include "G4TransportationProcessType.hh"
//TRANSPORTATION, COUPLED_TRANSPORTATION, STEP_LIMITER, USER_SPECIAL_CUTS, NEUTRON_KILLER
#include "G4StepStatus.hh"
//fWorldBoundary, fGeomBoundary, fAtRestDoItProc, fAlongStepDoItProc, fPostStepDoItProc, fUserDefinedLimit, fExclusivelyForcedProc, fUndefined
#include "G4TrackStatus.hh"
//fAlive, fStopButAlive, fStopAndKill, fKillTrackAndSecondaries, fSuspend, fPostponeToNextEvent

#include "MyDetectorConstruction.hh"
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wshadow"
#include "TH1F.h"
#pragma GCC diagnostic pop

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

//double MyG4BasedAnalysis::QE[14] = {45., 40., 35., 32., 29., 27., 22., 20., 14., 9., 4., 2., 0.1, 0.0};
MyG4BasedAnalysis::MyG4BasedAnalysis()
{
    SetFileName("g4output.root");

    //-------
    //#ANALYSIS 1. 初始化变量
    fTrkLen = 0;
    depositEnergyTot = 0;
    pTrans = 0;
    fEdeps.clear();
    fHitsX.clear();
    fHitsY.clear();
    fHitsZ.clear();
    eleDir = G4ThreeVector(.0, .0, .0);
    pHitsNum = 0;

    angleX = 0;
    MCtruthHit = 0;

    motherTrack = NULL;
    //double QE[14]={1,1,1,1,1,1,1,1,1,1,1,1,1,1};
    double QE[25] = {0.34461, 0.428418, 0.3675, 0.2852, 0.2212, 0.1756, 0.1360, 0.1025, 0.07218, 0.05549, 0.04489, 0.03431, 0.02372, 0.01618, 0.008638, 0.0056712, 0.0027077, 0.0028, 0.0013554, 0};
    //double QE[14] = {0.45, 0.40, 0.35, 0.32, 0.29, 0.27, 0.22, 0.20, 0.14, 0.09, 0.04, 0.02, 0.001, 0.0};
    fQE = new TH1F("fQE", "fQE", 25, 120, 210);
    for (int i = 0; i < 25; i++)
        fQE->SetBinContent(i, QE[i]);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

MyG4BasedAnalysis::~MyG4BasedAnalysis()
{
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void MyG4BasedAnalysis::BeginOfRunAction()
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::BeginOfRunAction()" << G4endl;

    auto analysisManager = G4AnalysisManager::Instance();
    // Default settings
    analysisManager->SetNtupleMerging(true);
    // Note: merging ntuples is available only with Root output

    analysisManager->SetVerboseLevel(1);
    analysisManager->OpenFile(fFileName);

    //-------
    //#ANALYSIS 2. 定义Ntuple结构

    // Creating 1D histograms
    //analysisManager->SetFirstHistoId(1);
    //analysisManager->CreateH1("phEng", "photon energy", 50, 0., 100); // h1 Id = 0

    // Creating 2D histograms
    //analysisManager->CreateH2("HitOnAnode", "Cherenkov photon hits on the anodes", // h2 Id = 0
    //                          50, -1000., 1000, 50, -1000., 1000.);

    // Creating ntuple
    //
    analysisManager->SetFirstNtupleId(1);
    analysisManager->CreateNtuple("CerenkovPosition", "Hits");
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("wavelen");
    analysisManager->FinishNtuple();

    //RICH，产生点，光波长、电子散射的分布

    analysisManager->CreateNtuple("boxPoint", "boxpoint");
    analysisManager->CreateNtupleDColumn("positionX");
    analysisManager->CreateNtupleDColumn("positionY");
    analysisManager->CreateNtupleDColumn("positionZ");
    analysisManager->CreateNtupleDColumn("wavelen2");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("photon", "Hits");
    analysisManager->CreateNtupleDColumn("number");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("hitposition", "hitposition");
    analysisManager->CreateNtupleDColumn("hitX");
    analysisManager->CreateNtupleDColumn("hitY");
    analysisManager->CreateNtupleDColumn("hitZ");
    analysisManager->CreateNtupleDColumn("wavelen3");
    analysisManager->FinishNtuple();
    /*
    analysisManager->CreateNtuple("ScatterTot","Hits");
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->FinishNtuple();
    */
    //QE waveLenth
    /*RICH beam test使用
    analysisManager->CreateNtuple("CerenkovLight", "Hits");

    analysisManager->CreateNtupleDColumn("X", fHitsX);
    analysisManager->CreateNtupleDColumn("Y", fHitsY);
    analysisManager->CreateNtupleDColumn("Z", fHitsZ);
    analysisManager->CreateNtupleDColumn("wavelen");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("photonNum", "Hits");
    analysisManager->CreateNtupleIColumn("pNum");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("Cerenkov","Hits");
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("MCtruth", "Hits");
    analysisManager->CreateNtupleDColumn("AngleX");
    analysisManager->CreateNtupleDColumn("AngleY");
    analysisManager->CreateNtupleDColumn("ANgleZ");
    analysisManager->FinishNtuple();
    */

    analysisManager->CreateNtuple("CerenkovLight", "Hits");
    analysisManager->CreateNtupleDColumn("X", fHitsX);
    analysisManager->CreateNtupleDColumn("Y", fHitsY);
    analysisManager->CreateNtupleDColumn("Z", fHitsZ);
    //analysisManager->CreateNtupleDColumn("wavelen");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("photonNum", "Hits");
    analysisManager->CreateNtupleIColumn("pNum");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("MCtruth", "Hits");
    analysisManager->CreateNtupleDColumn("angleX");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("Cerenkov", "Hits");
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("waveLen");
    analysisManager->FinishNtuple();
    /*准直的相关NTuple,version 1
    analysisManager->CreateNtuple("Gamma","Hits");
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("GammaEng");
    analysisManager->CreateNtupleDColumn("Px");
    analysisManager->CreateNtupleDColumn("Py");
    analysisManager->CreateNtupleDColumn("Pz");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("Electron","Hits");
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("EleEng");    
    analysisManager->CreateNtupleDColumn("Px");
    analysisManager->CreateNtupleDColumn("Py");
    analysisManager->CreateNtupleDColumn("Pz");
    analysisManager->CreateNtupleDColumn("parentID");
    analysisManager->CreateNtupleDColumn("PID");
    analysisManager->FinishNtuple(); 

    analysisManager->CreateNtuple("Neutron","Hits");
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");     
    analysisManager->CreateNtupleDColumn("NeuEng");   
    analysisManager->CreateNtupleDColumn("Px");
    analysisManager->CreateNtupleDColumn("Py");
    analysisManager->CreateNtupleDColumn("Pz");
    analysisManager->CreateNtupleDColumn("parentID");
    analysisManager->CreateNtupleDColumn("PID");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("GammaMid","Hits");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("GammaMidEng");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("EleMid","Hits");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("EleMidEng");    
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("NeuMid","Hits");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("NeuMidEng");    
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("GammaDis", "Hits");
    analysisManager->CreateNtupleDColumn("Gammatheta");
    analysisManager->CreateNtupleDColumn("Gammaphi");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("EleDis", "Hits");
    analysisManager->CreateNtupleDColumn("Eletheta");
    analysisManager->CreateNtupleDColumn("Elephi");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("NeuDis", "Hits");
    analysisManager->CreateNtupleDColumn("Neutheta");
    analysisManager->CreateNtupleDColumn("Neuphi");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("GammaCenter","Hits");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("GammaCenEng");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("EleCenter","Hits");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("EleCenEng");
    analysisManager->FinishNtuple();
    
    analysisManager->CreateNtuple("NeuCenter","Hits");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("NeuCenEng");
    analysisManager->FinishNtuple();
    */

    /*Collomator v2
    //10 cm 圆形
    analysisManager->CreateNtuple("Gamma","Hits");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("GammaEng");
    analysisManager->CreateNtupleDColumn("Px");
    analysisManager->CreateNtupleDColumn("Py");
    analysisManager->CreateNtupleDColumn("Pz");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("Electron","Hits");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("EleEng");
    analysisManager->CreateNtupleDColumn("Px");
    analysisManager->CreateNtupleDColumn("Py");
    analysisManager->CreateNtupleDColumn("Pz");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("Neutron","Hits");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("NeuEng");
    analysisManager->CreateNtupleDColumn("Px");
    analysisManager->CreateNtupleDColumn("Py");
    analysisManager->CreateNtupleDColumn("Pz");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("Particle","Hits");
    analysisManager->CreateNtupleSColumn("ParticleName");
    analysisManager->FinishNtuple();
    */

    /*
    analysisManager->CreateNtuple("Charged", "Hits"); // ntuple Id = 1
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("PX");
    analysisManager->CreateNtupleDColumn("PY");
    analysisManager->CreateNtupleDColumn("PZ");
    analysisManager->CreateNtupleDColumn("Process");
    analysisManager->CreateNtupleDColumn("ProcessSub");
    analysisManager->CreateNtupleDColumn("PID");
    analysisManager->CreateNtupleDColumn("ParentID");

    analysisManager->CreateNtuple("Neutral", "Hits"); // ntuple Id = 2
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("PX");
    analysisManager->CreateNtupleDColumn("PY");
    analysisManager->CreateNtupleDColumn("PZ");
    analysisManager->CreateNtupleDColumn("PID");
    analysisManager->CreateNtupleDColumn("ParentID");

    analysisManager->CreateNtuple("Track", "Hits"); // ntuple Id = 3
    analysisManager->CreateNtupleDColumn("TrkLen");
    analysisManager->CreateNtupleDColumn("Edeps", fEdeps);
    analysisManager->CreateNtupleDColumn("HitsX", fHitsX);
    analysisManager->CreateNtupleDColumn("HitsY", fHitsY);
    analysisManager->CreateNtupleDColumn("HitsZ", fHitsZ);

    analysisManager->CreateNtuple("NewTrk", "Hits"); // ntuple Id = 4
    analysisManager->CreateNtupleDColumn("Eng");
    */

    /*TPCsys保存的Ntuple信息
    analysisManager->CreateNtuple("Neutral", "Hits"); // ntuple Id = 1
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("ParentID");
    analysisManager->CreateNtupleDColumn("time");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("momentum", "Hits");
    analysisManager->CreateNtupleDColumn("PX");
    analysisManager->CreateNtupleDColumn("PY");
    analysisManager->CreateNtupleDColumn("PZ");
    analysisManager->CreateNtupleDColumn("PTot");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("angles", "Hits");
    analysisManager->CreateNtupleDColumn("theta");
    analysisManager->CreateNtupleDColumn("phi");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("Tot", "Hits");
    analysisManager->CreateNtupleDColumn("len");
    analysisManager->CreateNtupleDColumn("engDep");
    analysisManager->CreateNtupleDColumn("dedx");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("pad","Hits");
    analysisManager->CreateNtupleDColumn("ypad");
    analysisManager->CreateNtupleDColumn("zpad");    
    analysisManager->FinishNtuple();
    
    analysisManager->CreateNtuple("midTIME","Hits");
    analysisManager->CreateNtupleDColumn("Mt");
    analysisManager->FinishNtuple();
    */
    return;
}

void MyG4BasedAnalysis::EndOfRunAction()
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::EndOfRunAction()" << G4endl;

    //-------
    //#ANALYSIS 6. 在Run结束的时候将ntuple保存到文件
    auto analysisManager = G4AnalysisManager::Instance();
    analysisManager->Write();
    analysisManager->CloseFile();
    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void MyG4BasedAnalysis::BeginOfEventAction(const G4Event *)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::BeginOfEventAction()" << G4endl;

    //-------
    //#ANALYSIS 3. 初始化Event开始的参数
    // pTotMax = 0;
    // fTrkLen = 0;
    pHitsNum = 0;
    // depositEnergyTot = 0;
    // pTrans = 0;
    fEdeps.clear();
    fHitsX.clear();
    fHitsY.clear();
    fHitsZ.clear();
    eleDir = G4ThreeVector(.0, .0, .0);

    angleX = 0;
    MCtruthHit = 0;
    number = 0;

    motherTrack = NULL;

    return;
}

void MyG4BasedAnalysis::EndOfEventAction(const G4Event *)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::EndOfEventAction()" << G4endl;

    //-------
    //#ANALYSIS 5. 在Event结束的时候将数据保存到ntuple
    //G4cout<<"shashahahs"<<number<<G4endl;

    auto analysisManager = G4AnalysisManager::Instance();
    analysisManager->FillNtupleDColumn(3, 0, number);
    analysisManager->AddNtupleRow(3);

    //analysisManager->FillNtupleDColumn(3, 0, angleX / MCtruthHit);
    //analysisManager->AddNtupleRow(3);

    /*TPC保存径迹以及沉积能量和dE/dX
    if (fTrkLen != 0)
    {
        analysisManager->FillNtupleDColumn(4, 0, fTrkLen);
        analysisManager->FillNtupleDColumn(4, 1, depositEnergyTot);
        analysisManager->FillNtupleDColumn(4, 2, depositEnergyTot / fTrkLen);
        analysisManager->AddNtupleRow(4);
    }
    */

    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//此处是产生新的径迹
G4ClassificationOfNewTrack MyG4BasedAnalysis::ClassifyNewTrack(const G4Track *aTrack)
{
    if (!active)
        return fUrgent;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::ClassifyNewTrack()" << G4endl;

    // G4int parentID = aTrack->GetParentID();
    // const G4ParticleDefinition *particle = aTrack->GetParticleDefinition();
    // G4int pdgID = particle->GetPDGEncoding();
    // G4int charge = particle->GetPDGCharge();
    // if (charge == 0 && pdgID == 0 && parentID == 1){//QE
    //     G4String proName = aTrack->GetCreatorProcess()->GetProcessName();
    //     if (proName == "Cerenkov"){
    //         G4double waveLen = 0.0012398 / aTrack->GetKineticEnergy();
    //         if (waveLen > 215 || waveLen < 150)
    //         {
    //             return fKill;
    //         }
    //         double qe = fQE->Interpolate(waveLen);
    //         if (G4UniformRand() > qe)
    //         {
    //             return fKill;
    //         }
    //     }
    // }
    G4int parentID = aTrack->GetParentID();
    // G4Step *newStep = aTrack->GetStep();
    // auto *postStep = newStep->GetPostStepPoint();
    // auto *preStep = newStep->GetPreStepPoint();
    // if (postStep == NULL) return;
    // auto *postStepVolume = postStep->GetTouchableHandle()->GetVolume();
    // auto *preStepVolume = preStep->GetTouchableHandle()->GetVolume();
    // G4String nameNew = postVolume->GetName();

    // if(parentID == 0){
    //     // eleDir = aTrack->GetVertexMomentumDirection();
    //     //G4cout << "ele momDir is: (" << eleDir.x() << ", " << eleDir.y() << ", " << " )" << G4endl;
    // }
    const G4ParticleDefinition *particle = aTrack->GetParticleDefinition();
    G4int pdgID = particle->GetPDGEncoding();
    G4int charge = particle->GetPDGCharge();
    G4ThreeVector momDir = aTrack->GetMomentumDirection();
    //G4double eleX = 1;

    auto analysisManager = G4AnalysisManager::Instance(); //qe筛选以及计算MC truth
    if (charge == 0 && pdgID == 0 && parentID == 1)
    {
        eleDir = motherTrack->GetMomentumDirection();
        G4String proName = aTrack->GetCreatorProcess()->GetProcessName();
        if (proName == "Cerenkov")
        {
            G4double waveLen = 0.0012398 / aTrack->GetKineticEnergy();
            if (waveLen > 210 || waveLen < 120)
                return fKill;
            double qe = fQE->Interpolate(waveLen);
            if (G4UniformRand() > qe)
            {
                return fKill;
            }
            angleX += acos(momDir.x() * eleDir.x() + momDir.y() * eleDir.y() + momDir.z() * eleDir.z());
            MCtruthHit++;
            //G4cout << "MCtruthHit: " << MCtruthHit << G4endl;
            // G4cout << G4endl
            //        << "electron direction: " << "( " << eleDir.x() << " ," << eleDir.y() << " ,"
            //        << eleDir.z() << " )" << "norm is: " <<
            //        sqrt(eleDir.x()*eleDir.x()+eleDir.y()*eleDir.y()+eleDir.z()*eleDir.z())<< G4endl;
            // G4cout << "photon direction: " << "( " << momDir.x() << " ," << momDir.y() << " ,"
            //        << momDir.z() << " )" <<
            //        "norm of photon is: " << sqrt(momDir.x()*momDir.x()+momDir.y()*momDir.y()+momDir.z()*momDir.z())
            //        << G4endl;
            // G4ThreeVector vertex = aTrack->GetVertexMomentumDirection();
            // G4cout << "vertex info: " <<"( " << vertex.x() << ", " << vertex.y() << " ," << vertex.z() << " )" << G4endl;

            // analysisManager->FillNtupleDColumn(3,0,angleX);
            // analysisManager->AddNtupleRow(3);

            //G4cout<<"get track:LIUQ: "<<motherTrack->GetMomentumDirection()<<G4endl;
        }
    }
    //G4cout << G4endl;
    return fUrgent;
    //-------
    //#ANALYSIS 4.1 在生成新Track的时候保存相应数据

    //if (aTrack->GetParticleDefinition() == G4Gamma::Gamma())
    //    return fKill;

    //把QE放在前面，并且从粒子产生处获得MCtruth
    // G4int parentID = aTrack->GetParentID();
    // const G4ParticleDefinition *particle = aTrack->GetParticleDefinition();
    // G4int pdgID = particle->GetPDGEncoding();
    // G4int charge = particle->GetPDGCharge();
    // G4ThreeVector momDir = aTrack->GetMomentumDirection();
    // //G4double eleX = 1;
    // G4ThreeVector eleDir = G4ThreeVector(1.0, 0.0, 0.0);
    // if (parentID == 0){
    //     eleDir = momDir;
    // }
    // auto analysisManager = G4AnalysisManager::Instance();//qe筛选以及计算MC truth
    // if (charge == 0 && pdgID == 0 && parentID == 1){
    //     G4String proName = aTrack->GetCreatorProcess()->GetProcessName();
    //     if (proName == "Cerenkov"){
    //         G4double waveLen = 0.0012398 / aTrack->GetKineticEnergy();
    //         if (waveLen > 215 || waveLen < 150)
    //         {
    //             return fKill;
    //         }
    //         double qe = fQE->Interpolate(waveLen);
    //         if (G4UniformRand() > qe)
    //         {
    //             return fKill;
    //         }
    //         G4ThreeVector momDir = aTrack->GetMomentumDirection();
    //         G4double angleX = acos( momDir.x() * eleDir.x() + momDir.y() * eleDir.y()+momDir.z()*eleDir.z() );
    //         analysisManager->FillNtupleDColumn(3,0,angleX);
    //         analysisManager->AddNtupleRow(3);
    //     }
    // }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MyG4BasedAnalysis::PreTrackingAction(const G4Track *aTrack)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::PreTrackingAction()" << G4endl;
    //pN = 0;
    // G4ParticleDefinition *particle = aTrack->GetDefinition();
    // if (particle != NULL){
    // G4String name = particle->GetParticleName();
    // G4cout << name << G4endl;
    // if (name == "opticalphoton"){
    // G4String processName = aTrack->GetCreatorProcess()->GetProcessName();
    // G4cout <<  " processName: " << processName << G4endl;
    // }
    // }

    //通常用法为：
    //  保存某类粒子的相关信息(更多的时候是母粒子的次级粒子)，如产生点位置，产生点时间，物理过程，粒子种类等。
    //  特别是衰变的粒子或打靶实验
    /*
    {
        G4ParticleDefinition *particle = aTrack->GetDefinition();
        G4String name = particle->GetParticleName();
        G4int pid = particle->GetPDGEncoding();
        G4int Z = particle->GetAtomicNumber();
        G4int A = particle->GetAtomicMass();
        G4double charge = particle->GetPDGCharge();

        G4int ID = aTrack->GetTrackID();
        G4double energy = aTrack->GetKineticEnergy();
        G4double t0 = aTrack->GetGlobalTime();
        G4double weight = aTrack->GetWeight();
        G4int processType = aTrack->GetCreatorProcess()->GetProcessSubType(); //查看ProcessType的头文件定义
        G4String processName = aTrack->GetCreatorProcess()->GetProcessName();

        //还可根据坐标判断是否在Region内
        const G4ThreeVector& position = track->GetPosition();

        G4Region *region = G4RegionStore::GetInstance()->GetRegion("regName");
        int N =  region->GetNumberOfRootVolumes();
        std::vector<G4LogicalVolume*>::iterator it_logicalVolumeInRegion = region->GetRootLogicalVolumeIterator();

        bool inside_target = false;

        for(int i = 0; i < N ; i++, it_logicalVolumeInRegion++)
        {
            EInside test_status = (*it_logicalVolumeInRegion)->GetSolid()->Inside(position) ;
            if(test_status == kInside)
            {
                //...
                break;
            }
            else if (test_status == kSurface)
            {
                //...
                break;
            }
        }

        //...
    }
    */

    //-------
    //#ANALYSIS 4.2 在Tracking产生的时候保存相应数据

    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MyG4BasedAnalysis::PostTrackingAction(const G4Track *aTrack)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::PostTrackingAction()" << G4endl;
    // G4cout << pN << "<----pN in postTracking" << G4endl;
    // G4cout << G4endl;
    // G4cout << G4endl;

    //通常用法为：
    //  1. 根据粒子停止位置，计算能量泄漏等。
    //  2. 根据粒子停止时间，计算飞行时间等。
    //  3. 径迹总长度等信息
    /*
    {
        G4ParticleDefinition *particle = aTrack->GetDefinition();
        G4String name = particle->GetParticleName();
        G4int pid = particle->GetPDGEncoding();
        G4int Z = particle->GetAtomicNumber();
        G4int A = particle->GetAtomicMass();
        G4double charge = particle->GetPDGCharge();

        G4ThreeVector trkPos = aTrack->GetPosition();
        G4int ID = aTrack->GetTrackID();
        G4double energy = aTrack->GetKineticEnergy();
        G4double t1 = aTrack->GetGlobalTime();
        G4double weight = aTrack->GetWeight();
        G4int processType = aTrack->GetCreatorProcess()->GetProcessSubType(); //查看ProcessType的头文件定义
        G4String processName = aTrack->GetCreatorProcess()->GetProcessName();

        G4StepStatus stepStatus = aTrack->GetStep()->GetPostStepPoint()->GetStepStatus(); //检查最后一个step的状态, 常用fWorldBoundary判断leakage
        G4int nbSteps = aTrack->GetCurrentStepNumber();
        G4double Trleng = aTrack->GetTrackLength();

        G4TrackStatus trkStatus = aTrack->GetTrackStatus(); //检查track的状态，常用fStopAndKill判断
    }
    */

    //-------
    //#ANALYSIS 4.3 在Tracking终止的时候保存相应数据
    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MyG4BasedAnalysis::SteppingAction(const G4Step *aStep)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::SteppingAction()" << G4endl;

    //---
    //1. 相关参数的获取

    //1.1 Track的相关参数
    G4int parentID = aStep->GetTrack()->GetParentID();
    const G4ParticleDefinition *particle = aStep->GetTrack()->GetParticleDefinition();
    G4int pdgID = particle->GetPDGEncoding();
    G4int charge = particle->GetPDGCharge();
    // G4String particleName = particle->GetParticleName();

    if (parentID == 0 && motherTrack == NULL)
        motherTrack = aStep->GetTrack();

    //以下是G4Track常见的一些参数获取方法
    {
        /*
        G4int iStep = aTrack->GetCurrentStepNumber();
        G4int trackID = aTrack->GetTrackID();
        G4int parentID = aTrack->GetParentID();
        G4ParticleDefinition particle = aTrack->GetParticleDefinition();
        G4int pdgID = particle->GetPDGEncoding();
        G4String pdgName = particle->GetParticleName();
        G4cout << "-->" << iStep << ": tid=" << trackID << " " << pdgName << "(" << pdgID << ") parentID=" << parentID << G4endl;

        G4double energy = aTrack->GetKineticEnergy(); // kinetic energy at current step
        G4double engtot = aTrack->GetTotalEnergy();   // total energy (including m0)
        G4double velocity = aTrack->GetVelocity();    // mm/s
        G4ThreeVector moment = aTrack->GetMomentum();
        G4ThreeVector momDir = aTrack->GetMomentumDirection(); //unit vector
        G4ThreeVector polar = aTrack->GetPolarization();
        G4cout << "-->" << energy << " " << engtot << " " << velocity << " " << moment << " " << momDir << G4endl;

        //two different ways to get the process name
        G4String crtName = (parentID != 0) ? aTrack->GetCreatorModelName() : "none";
        G4String proName = (parentID != 0) ? aTrack->GetCreatorProcess()->GetProcessName() : "none";
        G4cout << "-->" << crtName << " " << proName << G4endl;

        G4ThreeVector trkPos = aTrack->GetPosition(); //current position of this track at the postStepPoint of this step
        G4double trkLen = aTrack->GetTrackLength(); //current track-length of this track at this postStepPoint of this step
        G4double gtime = aTrack->GetGlobalTime();
        G4double ltime = aTrack->GetLocalTime();
        G4double ptime = aTrack->GetProperTime();
        G4cout << "-->trkPos=" << trkPos << " " << trkLen << " " << gtime << " " << ltime << " " << ptime << G4endl;
        */
    }

    //1.2 Step的相关参数
    G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
    G4StepPoint *postStepPoint = aStep->GetPostStepPoint();

    //以下是G4Step常见的一些参数获取方法
    {
        /*
        G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
        G4StepPoint *postStepPoint = aStep->GetPostStepPoint();

        G4double stepLen = aStep->GetStepLength(); //the same as aTrack->GetTrackLength()
        G4double engDep = aStep->GetTotalEnergyDeposit(); // energy deposit inside this step
        G4double engDep2 = aStep->GetNonIonizingEnergyDeposit();
        G4bool isFirst = aStep->IsFirstStepInVolume();
        G4bool isLast = aStep->IsLastStepInVolume();
        G4cout << "-->" << stepLen << " " << engDep << " " << engDep2 << " " << isFirst << " " << isLast << G4endl;

        G4ThreeVector prePos = preStepPoint->GetPosition();
        G4double preGtime = preStepPoint->GetGlobalTime();
        G4double preLtime = preStepPoint->GetLocalTime();
        G4double prePtime = preStepPoint->GetProperTime();
        G4cout << "-->prePos=" << prePos << " " << preGtime << " " << preLtime << " " << prePtime << G4endl;
        G4ThreeVector postPos = postStepPoint->GetPosition(); // same as the information get from aTrack
        G4double postGtime = postStepPoint->GetGlobalTime();
        G4double postLtime = postStepPoint->GetLocalTime();
        G4double postPtime = postStepPoint->GetProperTime();
        G4cout << "-->postPos=" << postPos << " " << postGtime << " " << postLtime << " " << postPtime << G4endl;
        
        G4ThreeVector preMomDir = preStepPoint->GetMomentumDirection();
        G4ThreeVector preMomtum = preStepPoint->GetMomentum();
        G4double preEngTot = preStepPoint->GetTotalEnergy();
        G4double preEng = preStepPoint->GetKineticEnergy();
        G4double preV = preStepPoint->GetVelocity();
        G4double preBeta = preStepPoint->GetBeta();
        G4double preGamma = preStepPoint->GetGamma();

        G4ThreeVector postMomDir = postStepPoint->GetMomentumDirection();
        G4ThreeVector postMomtum = postStepPoint->GetMomentum();
        G4double postEngTot = postStepPoint->GetTotalEnergy();
        G4double postEng = postStepPoint->GetKineticEnergy();
        G4double postV = postStepPoint->GetVelocity();
        G4double postBeta = postStepPoint->GetBeta();
        G4double postGamma = postStepPoint->GetGamma();
        
        G4double mass = preStepPoint->GetMass();
        G4double charge = preStepPoint->GetCharge();
        G4double magmom = preStepPoint->GetMagneticMoment();

        auto *pVolume = postStepPoint->GetTouchableHandle()->GetVolume();
        auto *pLogicVolume = (pVolume!=NULL) ? pVolume->GetLogicalVolume() : NULL;

        */
    }

    //1.3 以下是拿DetectorConstruction相关参数的方法，根据具体情况修改即可
    {
        /*
        //拿Detector指针
        const MyDetectorConstruction* detectorConstruction = static_cast<const MyDetectorConstruction*> (G4RunManager::GetRunManager()->GetUserDetectorConstruction());

        //查找已知名字的逻辑体
        G4Box *fBox; 
        G4LogicalVolume* fLogicVol = G4LogicalVolumeStore::GetInstance()->GetVolume("YourLogicVolumeName");
        if ( fLogicVol ) fBox = dynamic_cast<G4Box*>(fLogicVol->GetSolid());
        fBox->GetXHalfLength() //这样就可以拿到具体的参数了。根据具体情况修改即可
        */
    }

    //-------
    //#ANALYSIS 4.4 在Steppinging的时候保存相应数据

    //---
    //2. 添加一些判断，并保存对应的数据。以下为演示，且按ANALYSIS 2. Ntuple定义的结构进行保存

    //G4String proName = postStepPoint->GetProcessDefinedStep()->GetProcessName();
    //G4ProcessType proType = postStepPoint->GetProcessDefinedStep()->GetProcessType();
    //G4int proSubType = postStepPoint->GetProcessDefinedStep()->GetProcessSubType();
    auto *preVolume = preStepPoint->GetTouchableHandle()->GetVolume();
    auto *postVolume = postStepPoint->GetTouchableHandle()->GetVolume();
    G4String PreVolume = preVolume->GetName();
    //为了获取前、后步所在的volume的名称
    if (postVolume == NULL)
        return;
    G4String name = particle->GetParticleName();
    G4String proName = postStepPoint->GetProcessDefinedStep()->GetProcessName();
    //G4int parentID = aTrack->GetParentID();
    G4LogicalVolume *preNowVolume = preVolume->GetLogicalVolume();
    G4LogicalVolume *postNowVolume = postVolume->GetLogicalVolume();
    G4ThreeVector postPos = postStepPoint->GetPosition();
    G4double optX = postPos.x();
    G4double optY = postPos.y();
    G4double optZ = postPos.z();
    G4ThreeVector postMom = postStepPoint->GetMomentum();
    G4double momX = postMom.x();
    G4double momY = postMom.y();
    G4double momZ = postMom.z();
    G4double momTot = sqrt(momX * momX + momY * momY + momZ * momZ);
    G4ThreeVector postMomDir = postStepPoint->GetMomentumDirection();
    G4ThreeVector preMomDir = preStepPoint->GetMomentumDirection();
    G4double preEng = preStepPoint->GetKineticEnergy();
    G4double postEng = postStepPoint->GetKineticEnergy();
    //G4double waveLen = 0.0012398 / aTrack->GetKineticEnergy();
    G4double waveLen = 0.001242375 / postEng;

    auto analysisManager = G4AnalysisManager::Instance();
    // if(postNowVolume->GetName() == "QuartzBoxVol" && preNowVolume->GetName() == "RICH" && parentID == 0){
    //     // G4cout << "this is the ele MomDir when into the quratz: (" << momX / momTot << ", " << momY / momTot << ", " << momZ / momTot << " )" << G4endl;
    //     // G4cout << "this is ele MomDir (" << postMomDir.x() << ", " << postMomDir.y() << ", " << postMomDir.z() << " )" << G4endl;
    //     eleDir = postMomDir;
    // }

    G4String particleName = aStep->GetTrack()->GetDynamicParticle()->GetDefinition()->GetParticleName();

    if (particleName == "opticalphoton") // for the 'opticalphoton'
    {
        //if (aStep->GetTrack()->GetTrackStatus() != fAlive)
        {

            G4ThreeVector start = aStep->GetTrack()->GetVertexPosition();
            G4ThreeVector momentum = aStep->GetTrack()->GetVertexMomentumDirection();
            G4ThreeVector final = aStep->GetTrack()->GetPosition();

            double x = final.x();
            double y = final.y();
            double z = final.z();
            double r = sqrt(x*x+y*y+z*z);
            if (fabs(x)>10)
            {
                G4cout<<"????"<<G4endl;
                analysisManager->FillNtupleDColumn(1, 0, final.x());
                analysisManager->FillNtupleDColumn(1, 1, final.y());
                analysisManager->FillNtupleDColumn(1, 2, final.z());
                analysisManager->FillNtupleDColumn(1, 3, waveLen);
                analysisManager->AddNtupleRow(1);

            }
            return;
        }
    }

    /*
    if (pdgID == 0 && postNowVolume->GetName() == "boardVol" && PreVolume == "ballPhys" && parentID == 1)
    {
        //if(pdgID==0 && postNowVolume->GetName() == "boardVol" && PreVolume=="ballPhys"&& parentID == 1){
        //G4double optEng = 0.0012398 / aTrack->GetKineticEnergy(); //convert to [nm]
        //for (int i = 0; i < 14; i++) fQE->SetBinContent(i + 1, QE[i]);
        //double qe = fQE->Interpolate(waveLen) / 100.;
        //if (G4UniformRand() > qe || (waveLen > 150 && waveLen < 220)){
        //G4String proName = aStep->GetCreatorProcess()->GetProcessName();
        //G4cout<<"shashasha"<<proName<<G4endl;
        //if (proName == "Cerenkov"){

        analysisManager->FillNtupleDColumn(1, 0, postPos.x());
        analysisManager->FillNtupleDColumn(1, 1, postPos.y());
        analysisManager->FillNtupleDColumn(1, 2, postPos.z());
        analysisManager->FillNtupleDColumn(1, 3, waveLen);
        analysisManager->AddNtupleRow(1);
        number = number + 1;
    }
    */


    //G4cout<<"hhhhhh"<<proName<<G4endl;
/*
    if (pdgID == 0 && parentID == 1)
    //if(pdgID==0 && parentID == 1)
    {
        if (preMomDir.x() * postMomDir.x() < 0 || preMomDir.y() * postMomDir.y() < 0 || preMomDir.z() * postMomDir.z() < 0)
        {
            //G4cout<<"hhh"<<G4endl;
            analysisManager->FillNtupleDColumn(2, 0, postPos.x());
            analysisManager->FillNtupleDColumn(2, 1, postPos.y());
            analysisManager->FillNtupleDColumn(2, 2, postPos.z());
            analysisManager->FillNtupleDColumn(2, 3, waveLen);
            analysisManager->AddNtupleRow(2);
        }
    }
*/

    // if (postNowVolume->GetName() == "FR4BoxVol" && pdgID == 0 && parentID == 1){
    //     fHitsX.push_back(optX);
    //     fHitsY.push_back((optY/cos(3.14159265359/6)+12.5));//hit point position = 8
    //     fHitsZ.push_back(optZ);
    //     pHitsNum++;
    //     //G4cout << "possition1 of " << name << "is ( " << postPos.x() << ", " << postPos.y() << ", " << postPos.z() << " )" << G4endl;
    //     //G4cout << "possition2 of " << name << "is ( " << optX << ", " << optY << ", " << optZ << " )" << G4endl;
    //     // analysisManager->FillNtupleDColumn(4,0,optX);
    //     // analysisManager->FillNtupleDColumn(4,1,(optY/cos(3.1415926/6)+12.5)/5);
    //     // analysisManager->FillNtupleDColumn(4,2,optZ/5+16);
    //     // analysisManager->FillNtupleDColumn(4,3,waveLen);
    //     // analysisManager->AddNtupleRow(4);
    // }

    // fHitsX.clear();
    // fHitsY.clear();//hit point position = 8
    // fHitsZ.clear();
    /*这个Block是切伦科夫光的相关信息
    if (postNowVolume->GetName() == "QuartzBoxVol" && proName == "Cerenkov" && parentID == 0){
        analysisManager->FillNtupleDColumn(1,0,optX);
        analysisManager->FillNtupleDColumn(1,1,optY);
        analysisManager->FillNtupleDColumn(1,2,optZ);
        analysisManager->AddNtupleRow(1);
    }
    if (preNowVolume->GetName() == "QuartzBoxVol" && postNowVolume->GetName() == "ArBoxVol" && parentID == 0){
        analysisManager->FillNtupleDColumn(3,0,optX);
        analysisManager->FillNtupleDColumn(3,1,optY);
        analysisManager->FillNtupleDColumn(3,2,optZ);
        analysisManager->AddNtupleRow(3);
    }
    if (postNowVolume->GetName() == "FR4boxVol" && parentID == 0){
        analysisManager->FillNtupleDColumn(4,0,optX);
        analysisManager->FillNtupleDColumn(4,1,optY);
        analysisManager->FillNtupleDColumn(4,2,optZ);
        analysisManager->AddNtupleRow(4);
    }
    */
    //QE considered
    // if (pdgID == 0){
    // pN++;
    // }

    // if (postNowVolume->GetName() == "FR4BoxVol" && pdgID == 0)
    // {
    //    if (waveLen > 215 || waveLen < 150)
    //     {
    //         return;
    //     }

    //     double qe = fQE->Interpolate(waveLen);
    //     if (G4UniformRand() > qe)
    //         return;
    //     fHitsX.push_back(optX);
    //     fHitsY.push_back(optY/cos(3.14159265359/6));//hit point position = 8
    //     fHitsZ.push_back(optZ);
    //     analysisManager->FillNtupleDColumn(1, 3, waveLen);
    //     analysisManager->AddNtupleRow(1);
    //     // G4cout << particleName << G4endl;
    //     pHitsNum++;
    //     //G4cout << pHitsNum << "<----pHitsNum " << G4endl;
    // }
    // analysisManager->FillNtupleDColumn(3,0,optX);
    // analysisManager->FillNtupleDColumn(3,1,(optY/cos(3.14159265359/6)+5)/5);
    // analysisManager->FillNtupleDColumn(3,2,optZ/5+16);
    // analysisManager->AddNtupleRow(3);

    // }
    // if (pdgID == 0 && iStep == 0 && postNowVolume->GetName() == "QuartzBoxVol" && parentID == 1){
    // G4String processName = aTrack->GetCreatorProcess()->GetProcessName();
    // double angleX = 0;
    // double angleY = 0;
    // double angleZ = 0;
    // if (processName == "Cerenkov"){
    // angleX = acos( postMomDir.x() );
    // angleY = acos( postMomDir.y() );
    // angleZ = acos( postMomDir.z() );
    // analysisManager->FillNtupleDColumn(4,0,angleX);
    // analysisManager->FillNtupleDColumn(4,1,angleY);
    // analysisManager->FillNtupleDColumn(4,2,angleZ);
    // analysisManager->AddNtupleRow(4);
    // }
    //G4cout << iStep << "<-------iStep of opticalphoton" << G4endl;
    //G4cout << G4endl;
    // }

    // if (pdgID == 0){
    // G4String processName = aTrack->GetCreatorProcess()->GetProcessName();
    // G4cout << processName << G4endl;
    // G4cout << "pN in Step: " << pN << G4endl;
    // G4cout << "step ID: " << iStep << G4endl;
    // G4cout << preNowVolume->GetName() <<G4endl;
    // G4cout << postNowVolume->GetName() << G4endl;
    // G4cout << G4endl;
    // }
    //
    // G4cout << "Name is :----->" << name << G4endl;
    // G4cout << "pdgID is:------>" << pdgID << G4endl;
    // G4cout << "ProName is:---->" << proName << G4endl;
    // G4cout << "Present Volume is:------>" << postNowVolume->GetName() << G4endl;
    // G4cout << "X position is:--->" << optX << G4endl;
    // G4cout << "parentID is --->" << parentID << G4endl;
    // G4cout << "IsFirstStepInVolume---->" << isFirst << G4endl;
    // G4cout << G4endl;

    /*准直version 1
    if (pdgID == 22 ){
        if(postNowVolume->GetName() == "ReadPadVol"){
            analysisManager->FillNtupleDColumn(1, 0, optX);
            analysisManager->FillNtupleDColumn(1, 1, optY);
            analysisManager->FillNtupleDColumn(1, 2, optZ);
            analysisManager->FillNtupleDColumn(1, 3, preEng);
            analysisManager->FillNtupleDColumn(1, 4, momX);
            analysisManager->FillNtupleDColumn(1, 5, momY);
            analysisManager->FillNtupleDColumn(1, 6, momZ);
            analysisManager->FillNtupleDColumn(1, 7, parentID);
            analysisManager->FillNtupleDColumn(1, 8, pdgID);
            analysisManager->AddNtupleRow(1);        
            if (optY > -50 && optY < 50 && optZ > -50 && optZ < 50){
                analysisManager->FillNtupleDColumn(4,0,optY);
                analysisManager->FillNtupleDColumn(4,1,optZ);
                analysisManager->FillNtupleDColumn(4,2,preEng);
                analysisManager->AddNtupleRow(4);
                double theta = acos(momX / momTot);
                double phi = acos((momY / momTot) / sin(theta));
                analysisManager->FillNtupleDColumn(7,0,theta);
                analysisManager->FillNtupleDColumn(7,1,phi);
                analysisManager->AddNtupleRow(7);               
            }
            if (optY > -15 && optY < 15 && optZ > -15 && optZ < 15){
                analysisManager->FillNtupleDColumn(10,0,optY);
                analysisManager->FillNtupleDColumn(10,1,optZ);
                analysisManager->FillNtupleDColumn(10,2,preEng);
                analysisManager->AddNtupleRow(10);
            }

        }
    }
    if (pdgID == 11 ){
        if (postNowVolume->GetName() == "ReadPadVol"){
            analysisManager->FillNtupleDColumn(2, 0, optX);
            analysisManager->FillNtupleDColumn(2, 1, optY);
            analysisManager->FillNtupleDColumn(2, 2, optZ);
            analysisManager->FillNtupleDColumn(2, 3, preEng);
            analysisManager->FillNtupleDColumn(2, 4, momX);
            analysisManager->FillNtupleDColumn(2, 5, momY);
            analysisManager->FillNtupleDColumn(2, 6, momZ);
            analysisManager->FillNtupleDColumn(2, 7, parentID);
            analysisManager->FillNtupleDColumn(2, 8, pdgID);
            analysisManager->AddNtupleRow(2);
            if (optY > -50 && optY < 50 && optZ > -50 && optZ < 50){
                analysisManager->FillNtupleDColumn(5,0,optY);
                analysisManager->FillNtupleDColumn(5,1,optZ);
                analysisManager->FillNtupleDColumn(5,2,preEng);
                analysisManager->AddNtupleRow(5);
                double theta = acos(momX / momTot);
                double phi = acos((momY / momTot) / sin(theta));
                analysisManager->FillNtupleDColumn(8,0,theta);
                analysisManager->FillNtupleDColumn(8,1,phi);
                analysisManager->AddNtupleRow(8);               
            }
            if (optY > -15 && optY < 15 && optZ > -15 && optZ < 15){
                analysisManager->FillNtupleDColumn(11,0,optY);
                analysisManager->FillNtupleDColumn(11,1,optZ);
                analysisManager->FillNtupleDColumn(11,2,preEng);
                analysisManager->AddNtupleRow(11);
            }
        }

    }
    if (pdgID == 2112 ){
        if (postNowVolume->GetName() == "ReadPadVol"){
            analysisManager->FillNtupleDColumn(3, 0, optX);
            analysisManager->FillNtupleDColumn(3, 1, optY);
            analysisManager->FillNtupleDColumn(3, 2, optZ);
            analysisManager->FillNtupleDColumn(3, 3, preEng);
            analysisManager->FillNtupleDColumn(3, 4, momX);
            analysisManager->FillNtupleDColumn(3, 5, momY);
            analysisManager->FillNtupleDColumn(3, 6, momZ);
            analysisManager->FillNtupleDColumn(3, 7, parentID);
            analysisManager->FillNtupleDColumn(3, 8, pdgID);
            analysisManager->AddNtupleRow(3);
            if (optY > -50 && optY < 50 && optZ > -50 && optZ < 50){
                analysisManager->FillNtupleDColumn(6,0,optY);
                analysisManager->FillNtupleDColumn(6,1,optZ);
                analysisManager->FillNtupleDColumn(6,2,preEng);
                analysisManager->AddNtupleRow(6);
                double theta = acos(momX / momTot);
                double phi = acos((momY / momTot) / sin(theta));
                analysisManager->FillNtupleDColumn(9,0,theta);
                analysisManager->FillNtupleDColumn(9,1,phi);
                analysisManager->AddNtupleRow(9);          
            }
            if (optY > -15 && optY < 15 && optZ > -15 && optZ < 15){
                analysisManager->FillNtupleDColumn(12,0,optY);
                analysisManager->FillNtupleDColumn(12,1,optZ);
                analysisManager->FillNtupleDColumn(12,2,preEng);
                analysisManager->AddNtupleRow(12);
            }
        }
    }*/

    /*准直器version 2
    if( pdgID == 22 && postNowVolume->GetName() == "ReadPadVol"){
        analysisManager->FillNtupleDColumn(1, 0, optY);
        analysisManager->FillNtupleDColumn(1, 1, optZ);
        analysisManager->FillNtupleDColumn(1, 2, preEng);
        analysisManager->FillNtupleDColumn(1, 3, momX);
        analysisManager->FillNtupleDColumn(1, 4, momY);
        analysisManager->FillNtupleDColumn(1, 5, momZ);
        analysisManager->AddNtupleRow(1);
    }
    if( pdgID == 11 && postNowVolume->GetName() == "ReadPadVol"){
        analysisManager->FillNtupleDColumn(2, 0, optY);
        analysisManager->FillNtupleDColumn(2, 1, optZ);
        analysisManager->FillNtupleDColumn(2, 2, preEng);
        analysisManager->FillNtupleDColumn(2, 3, momX);
        analysisManager->FillNtupleDColumn(2, 4, momY);
        analysisManager->FillNtupleDColumn(2, 5, momZ);
        analysisManager->AddNtupleRow(2);
    }            
    if( pdgID == 2112 && postNowVolume->GetName() == "ReadPadVol"){
        analysisManager->FillNtupleDColumn(3, 0, optY);
        analysisManager->FillNtupleDColumn(3, 1, optZ);
        analysisManager->FillNtupleDColumn(3, 2, preEng);
        analysisManager->FillNtupleDColumn(3, 3, momX);
        analysisManager->FillNtupleDColumn(3, 4, momY);
        analysisManager->FillNtupleDColumn(3, 5, momZ);
        analysisManager->AddNtupleRow(3);
    }
    if( postNowVolume->GetName() == "ReadPadVol"){
        analysisManager->FillNtupleSColumn(4,0,name);
        analysisManager->AddNtupleRow(4);
    }
    */

    /*
    if( preNowVolume->GetName() == "CollomatorVol" && postNowVolume->GetName() == "World"){
        analysisManager->FillNtupleSColumn(4,0,name);
        analysisManager->AddNtupleRow(4);
    }
    */

    /*这个block里面包含的是TPC所要保存的项目,共5个Ntuple，分别保存径迹、出射角分布、进入TPC粒子动量、坐标变换后投影坐标和中心时间
        if (parentID != 0){
            return;
        }
        if (preNowVolume->GetName() == "AluminumSERQianVol" && postNowVolume->GetName() == "TPCsys")
        {
            double theta = acos(momX / momTot);
            double phi = acos((momY / momTot) / sin(theta));
            analysisManager->FillNtupleDColumn(3, 0, theta);
            analysisManager->FillNtupleDColumn(3, 1, phi);
            analysisManager->AddNtupleRow(3);
        }
        if (preNowVolume->GetName() == "ArBoxTPCVol" && postNowVolume->GetName() == "ArBoxTPCVol")
        {
            analysisManager->FillNtupleDColumn(1, 0, optX);
            analysisManager->FillNtupleDColumn(1, 1, optY);
            analysisManager->FillNtupleDColumn(1, 2, optZ);
            analysisManager->FillNtupleDColumn(1, 3, parentID);
            G4double postGtime = postStepPoint->GetGlobalTime();
            analysisManager->FillNtupleDColumn(1, 4, postGtime);
            double engDep = aStep->GetTotalEnergyDeposit();
            double trkLen = aStep->GetStepLength();
            fTrkLen += trkLen;
            depositEnergyTot += engDep;
            //G4cout<<"---->"<<depositEnergyTot<<G4endl;
            analysisManager->AddNtupleRow(1);
            //G4cout << optX << "<-------" << G4endl;
            optX += -70+37.3982;
            optY += 87.7803;
            double yPad = sin(0.593)*optX + cos(0.593)*optY;
            analysisManager->FillNtupleDColumn(5,0,yPad);
            analysisManager->FillNtupleDColumn(5,1,optZ);
            analysisManager->AddNtupleRow(5);
            if (yPad > 75 && yPad < 85 && optZ > -5 && optZ < 5){
                analysisManager->FillNtupleDColumn(6,0,postGtime);
                analysisManager->AddNtupleRow(6);
            }
            

            
        }
        if (preNowVolume->GetName() == "TPCVol" && postNowVolume->GetName() == "ArBoxTPCVol")
        {
            analysisManager->FillNtupleDColumn(2, 0, momX);
            analysisManager->FillNtupleDColumn(2, 1, momY);
            analysisManager->FillNtupleDColumn(2, 2, momZ);
            analysisManager->FillNtupleDColumn(2, 3, momTot);
            analysisManager->AddNtupleRow(2);
        }

            
        G4double stepLen = aStep->GetStepLength();
        G4double engDep = aStep->GetTotalEnergyDeposit();
        fTrkLen += stepLen;
        depositEnergyTot += engDep;
        */

    /*fEdeps.push_back(aStep->GetTotalEnergyDeposit());
        fHitsX.push_back(postPos.x());
        fHitsY.push_back(postPos.y());
        fHitsZ.push_back(postPos.z());

        if (presentVolume->GetName() == "FR4BoxVol"){
        G4ThreeVector postPos = postStepPoint->GetPosition();
        G4double optX = postPos.x();
        G4double optY = postPos.y();
        G4double optZ = postPos.z();
        G4ThreeVector postMom = postStepPoint->GetMomentum();
        G4double momX = postMom.x();
        G4double momY = postMom.y();
        G4double momZ = postMom.z();
        G4double momTot = sqrt(momX*momX + momY*momY + momZ*momZ);
        auto analysisManager = G4AnalysisManager::Instance();  
        analysisManager->FillNtupleDColumn(2, 0, optX);
        analysisManager->FillNtupleDColumn(2, 1, optY);
        analysisManager->FillNtupleDColumn(2, 2, optZ);
        analysisManager->FillNtupleDColumn(2, 3, momX);
        analysisManager->FillNtupleDColumn(2, 4, momY);
        analysisManager->FillNtupleDColumn(2, 5, momZ);
        analysisManager->FillNtupleDColumn(2, 6, momTot);
        analysisManager->FillNtupleDColumn(2, 7, parentID);
        analysisManager->AddNtupleRow(2); }        
        */
}
/*
    //Ntuple1: 保存次级粒子中带电粒子信息
    if (parentID != 0 && charge != 0)
    {
        //if (aTrack->GetTrackStatus() != fStopAndKill) //只要track停止时的信息
        //    return;

        //G4cout << "==>"<<proName << G4endl;
        //if (proName != "eIoni") //只要电离过程
        //    return;

        auto analysisManager = G4AnalysisManager::Instance();
        analysisManager->FillNtupleDColumn(1, 0, postPos.x());
        analysisManager->FillNtupleDColumn(1, 1, postPos.y());
        analysisManager->FillNtupleDColumn(1, 2, postPos.z());
        analysisManager->FillNtupleDColumn(1, 3, postMomtum.x());
        analysisManager->FillNtupleDColumn(1, 4, postMomtum.y());
        analysisManager->FillNtupleDColumn(1, 5, postMomtum.z());
        analysisManager->FillNtupleDColumn(1, 6, double(proType));
        analysisManager->FillNtupleDColumn(1, 7, double(proSubType));
        analysisManager->FillNtupleDColumn(1, 8, double(pdgID));
        analysisManager->FillNtupleDColumn(1, 9, double(parentID));
        analysisManager->AddNtupleRow(1);
    }

    //Ntuple2: 保存中性粒子信息：
    if (charge == 0 && (parentID == 1 && pdgID == 0)) //要求来自入射粒子，且是光子
    {
        if (proName != "Cerenkov") //只要切伦科夫过程
            return;

        auto *pVolume = postStepPoint->GetTouchableHandle()->GetVolume();
        if (pVolume == NULL)
            return;

        G4LogicalVolume *presentVolume = pVolume->GetLogicalVolume();
        if (presentVolume->GetName() != "FR4BoxVol") //只要光子打到阳极板上的情况
            return;

        G4double optEng = 0.0012398 / aTrack->GetKineticEnergy(); //convert to [nm]
        G4double optX = postPos.x();
        G4double optY = postPos.y();
        G4double optZ = postPos.z();

        auto analysisManager = G4AnalysisManager::Instance();
        analysisManager->FillNtupleDColumn(0, 0, optEng);
        analysisManager->FillNtupleDColumn(0, 1, optX);
        analysisManager->FillNtupleDColumn(0, 2, optY);
        analysisManager->FillNtupleDColumn(0, 3, optZ);
        analysisManager->FillNtupleDColumn(0, 4, aTrack->GetVertexPosition().x());
        analysisManager->FillNtupleDColumn(0, 5, aTrack->GetVertexPosition().y());
        analysisManager->FillNtupleDColumn(0, 6, aTrack->GetVertexPosition().z());
        analysisManager->AddNtupleRow(2);
    }

    //Ntuple3: a)TrkLen在这里累加后，在EndOfEvent处保存； b)用vector来保存保存入射粒子的hit信息
    if (parentID == 0)
    {
        G4double stepLen = aStep->GetStepLength();
        fTrkLen += stepLen;

        fEdeps.push_back(aStep->GetTotalEnergyDeposit());
        fHitsX.push_back(postPos.x());
        fHitsY.push_back(postPos.y());
        fHitsZ.push_back(postPos.z());
    }

    return;
       if (parentID == 0)
    {
        auto *pVolume = postStepPoint->GetTouchableHandle()->GetVolume();
        if (pVolume == NULL)
            return;
        
        G4LogicalVolume *presentVolume = pVolume->GetLogicalVolume();
        if (presentVolume->GetName() == "ArBoxTPCVol")
        {
        G4ThreeVector postPos = postStepPoint->GetPosition();
        G4double optX = postPos.x();
        G4double optY = postPos.y();
        G4double optZ = postPos.z();

        G4ThreeVector postMom = postStepPoint->GetMomentum();
        G4double momX = postMom.x();
        G4double momY = postMom.y();
        G4double momZ = postMom.z();

        fEdeps.push_back(aStep->GetTotalEnergyDeposit());
        fHitsX.push_back(optX);
        fHitsY.push_back(optY);
        fHitsZ.push_back(optZ);
        G4double stepLen = aStep->GetStepLength();
        fTrkLen += stepLen;
        //G4cout << "the momentum is: " << postMom << G4endl;        
        auto analysisManager = G4AnalysisManager::Instance();        
        analysisManager->FillNtupleDColumn(1, 0, optX);
        analysisManager->FillNtupleDColumn(1, 1, optY);
        analysisManager->FillNtupleDColumn(1, 2, optZ);
        analysisManager->FillNtupleDColumn(1, 3, momX);
        analysisManager->FillNtupleDColumn(1, 4, momY);
        analysisManager->FillNtupleDColumn(1, 5, momZ);
        analysisManager->FillNtupleDColumn(1, 6, parentID);
        analysisManager->AddNtupleRow(1);
        }
    }
    */
